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Abstract 

The response of ecosystems to pertm'bations is considered from a thermo- 
dynamic perspective by acknowledging that, as for aU macroscopic systems and 
processes, the dynamics and stabihty of ecosystems is subject to definite ther- 
modynamic law. For open ecosystems, exchanging energy, work, and mass with 
the environment, the thermodynamic criteria come from non-equilibrium or irre- 
versible thermodynamics. For ecosystems during periods in which the boundary 
conditions may be considered as being constant, it is shown that criteria from 
irreversible thermodynamic theory are sufficient to permit a quantitative predic- 
tion of ecosystem response to perturbation. This framework is shown to provide 
a new perspective on the population dynamics of real ecosystems. The formal- 
ism is applied to the problem of the population oscillations of the southern pine 
beetle. 
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I. Introduction 

Most ecosystems are under considerable stress, having been perturbed by 
human intrusions including, population reduction or complete species annihila- 
tion, the introduction of foreign species, habitat destruction or fragmentation, 
contamination, and general global warming. We are now witnessing one of the 
highest rates of species disappearance and ecosystem collapse over the entire 
history of life on Earth. 

Fortunately, ecosystems can often recover from perturbations and can even 
evolve and adapt to new boundary conditions (Norberg et al., 2001). However, 
successful recovery depends on the inherent stability of the system, which is 
a complex function of the individual interactions among all the participating 
species and among species and their environment. Given that typical ecosys- 
tems contain over 3000 species (Polls, 1991), understanding the nature of this 
stability, and thus predicting ecosystem response to perturbation, is far from 
trivial, but indispensable for obtaining a quantitative understanding of ecosys- 
tem dynamics. 

Predicting ecosystem response to perturbation is, therefore, one of the most 
scientifically taxing yet important questions of our time. Present ecosystem dy- 
namics theory is based on empirically inspired but essentially ad hoc equations 
incorporating one-body parameters, such as inherent birth and death rates, and 
two-body effects of one species population on another through coefficients rep- 
resenting competition, predator-prey, symbiosis, or neutral interaction, as well 
as species-environment effects, incorporated through parameters such as the 
"carrying capacity" . This two-body "community matrix" approach, although 
widely recognized for its usefulness in revealing the general spectrum of the 
dynamics of model ecosystems (May, 1974), has had little application to pre- 
dicting real ecosystem response to perturbation. This is primarily due to the 
fact that the dynamical equations are ad hoc and that the community matrix is 
obtained from fits to time-series population data, and therefore can be expected 
to be representative of nature only within the limited range of the available 
population data. 

A further problem debilitating the community matrix approach is that it 
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is a two-body approach while species interactions are really of a many-body 
natnrc. These many-body effects are usually absorbed within so called "en- 
vironmental factors" which arc included in the dynamical equations as fitted 
constants. However, these "constants" are not really constant for perturbed 
ecosystems and anyhow fail to endow the resultant 2-body equations with the 
true dynamics of inherently many-body natural ecosystems. 

A clear indication of the failure of traditional theory is the fact that today's 
ecosystem health is usually surmised by making painstaking field counts of the 
populations of particular species and then, rather arbitrarily, deciding whether 
or not to include those species on an "in danger of extinction" list. Such one 
dimensional and last minute vigilance of ecosystems is not satisfactory for a 
number of reasons: first, it fails to treat the ecosystem as an integrated whole 
and could thereby miscalculate the gravity of the situation about to unfold; 
second, since many ecosystems have a natural cyclical, or even chaotic, but 
stable dynamics, it may be difficult to distinguish normal, but stable, periodic 
or chaotic behavior from a dangerous fall toward extinction; third, our human 
perspective tends to focus on species in which the individuals are physically 
large, easily observable, or likable, but not necessarily those key species that are 
most important to the stability of an ecosystem. Most important, however, is 
the fact that present ecosystem theory provides no information for designing an 
integral solution for arresting an impending catastrophe, other than, perhaps, 
suggesting that the endangered specie be protected by law. 

There is clearly a need for a more quantitative approach to population mod- 
eling based on fundamental science and measurements that can predict ecosys- 
tem dynamics for regions in population space for which no data exists. There 
has been a growing realization that such a quantitative theory of ecosystems 
will have a thermodynamic basis (Odum, 1969, Gallucci, 1973, Swenson, 1989, 
Michaelian, 2005). The reasons are compelling: First, thermodynamic laws 
derive from symmetry principles inherent in nature and thus are universal, ap- 
plicable in suitable form to all macroscopic systems and processes, irrespective 
of the types of interactions involved. Second, thermodynamics deals with a 
much reduced set of macroscopic variables which can be related with measur- 
able ecosystem variables (e.g. populations) involved in the dynamical patterns 
observed in Nature. Third, a number of stubborn problems and paradoxes ex- 
isting in traditional ecosystem theory appear to have a simple resolution in 



3 



terms of thermodynamic directives (Swenson, 1989, Michaelian, 2005). The 

objective of this paper is to demonstrate that, for ecosystems under constant 
boundary conditions, a non-equihbrium thermodynamic framework for the pop- 
ulation dynamics can lead to explicit predictions concerning ecosystem response 
to perturbation. 

In the following section we briefly outline the thermodynamic framework 
for treating ecosystems which has been presented elsewhere (Michaelian, 2005). 
In section 3 we present a simple model ecosystem and demonstrate how its 
population dynamics and stability characteristics are determined by thermody- 
namic constraints and criteria relating to energy, work, and mass flow among 
the species populations, and with the external environment. In section 4 we 
perturb this ecosystem and analyze the response as predicated on the basis of 
non-equilibrium thermodynamic formalism. Finally, in section 5, we discuss how 
our thermodynamic framework may have relevance in explaining the particular 
population dynamics observed in a real ecosystem; the outbreak of sustained 
population oscillations of the southern pine beetle. 

II. Thermodynamic Framework 

To avoid misinterpretation, it is prudent to first make a clear distinction 
between two existing but fundamentally different thermodynamic frameworks. 
Equilibrium thermodynamics deals with isolated systems and the fundamental 
state variable governing the evolution of the isolated system toward the stable 
equilibrium state is the total entropy, S. Irreversible thermodynamics deals 
with open systems or processes, such as ecosystems, which exchange energy, 
work, and mass between component parts and with the environment. Here, the 
variable governing the evolution toward the stable stationary state (for constant 
boundary conditions) is the time variation of the total entropy of the system, 
dS/ dt. Our framework is based on the latter irreversible thermodynamics and 
we employ only that part of this framework, known as classical, developed by 
Lars Onsager (Onsager, 1931) and lUya Prigogine (Prigogine, 1967) that has 
been verified empirically. 

As for any open system, the time variation of the total entropy of the ecosys- 
tem may be divided into a part due to the internal entropy production arising 
from irreversible processes occurring within the ecosystem itself, and a second 
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part due to the flow of entropy into, or out of, the ecosystem from the external 
environment (Prigogine, 1967), 

dS diS deS 

— = — h (1) 

dt dt dt ^ ' 

All macroscopic systems and processes, including ecosystems, are subject to 
definite thermodynamic law. The primary among these is the second law of 
thermodynamics which states that the internal production of entropy due to 
irreversible processes occurring within the system must be positive definite, 

For the case of ecosystems under the condition of constant external con- 
straints (see (Michaelian, 2005) for justification of this condition for a large 
class of ecosystems) classical irreversible thermodynamic theory states (Pri- 
gogine, 1967) that the system will eventually arrive at a thermodynamic sta- 
tionary state in which all macroscopic variables, including the total entropy, are 
stationary in time, 

Therefore, from ([T]), at the stationary state, 
implying from Eq. ^ that 

Maintaining an ecosystem in a stable thermodynamic stationary state thus re- 
quires a continuous negative flow of entropy into the system. This has already 
been emphasized by Schrodinger (Schrodinger, 1944), but was first recognized 
by Boltzmann (Boltzmann, 1886). 

The internal entropy production diS/dt can be written as a sum of general- 
ized thermodynamic forces X multiplied by generalized thermodynamic flows J 
(Prigogine, 1967) (for example, Xq = V (^)gradient of the inverse temperature, 
and Jq = heat flow), 

a 

The separation of the entropy production into its components of thermodynamic 
forces and flows is somewhat arbitrary and can often be chosen for convenience 
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in resolving the particular problem at hand. However, there are a number 
of conditions that must be met for any particular choice. The first condition 
is that the product of the force and flow gives units of entropy production, 
and the second is that symmetry aspects must be respected, for example, a 
scaler force cannot give rise to a vector flow (Katchalsky and Curran, 1975). 
We have shown (Michaelian, 2005) that ecosystem dynamics can be treated 
consistently within this irreversible thermodynamic framework by assigning the 
generalized thermodynamic forces to the species populations {Xa = Pa) (where 
a represents the species type) and the generalized flows to the flows of entropy 
( Jq, = diSa/dt) (due to flows of energy, work, or mass, to or from species a, see 
below) . 

Within this framework, it was shown (Michaelian, 2005) that the ecological 
steady state, prevalent in nature (Goldwasser & Roghgarden, f993; Polls, 1991), 
has the stability characteristics of the thermodynamic stationary state. In view 
of this, we have made the formal assertion that the ecological steady state 
is just a particular case of the more general thermodynamic stationary state 
(Michaehan, 2005). 

A further criterion from classical irreversible thermodynamic theory, con- 
sidered by Prigogine as the most general result of irreversible thermodynamic 
theory, and valid for constant external constraints, is that the rate of change 
of the internal entropy production, due to changes in the generalized forces X 
(the populations), is negative semi-definite; the general evolutionary criterium 
(Prigogine 1967), 

— ; — < where V = —— . (7) 

dt dt 

Equation ([7]) implies that, under constant boundary conditions, all natural 
changes in the species populations must be in such a manner so as to reduce 
the internal production of entropy. This is a powerful auxiliary criterion on 
ecosystem response to perturbation and it will be shown below that this, to- 
gether with the second law of thermodynamics, and the fact that a system with 
constant external constraints must arrive at a thermodynamic stationary state, 
effectively determines the population dynamics that the ecosystem can assume. 
In this way, we can predict the actual dynamic response of the ecosystem to 
perturbation, be it either toward recovery, toward a new dynamics, or toward 
extinction. 



6 



III. Model Ecosystem 



We now present our thermodynamic framework for an over simplified but 
illustrative 3-species model ecosystem, including up to 3-body interaction terms. 
Two of the populations, pi and P2, are considered variable, while the third, p^, 
is fixed, and represents the constant boundary conditions over the ecosystem; 
such as the constant supply of nutrients due to a primary producer species. For 
example, in our application to the pine beetle ecosystem, to be detailed in the 
section 5, p3 would represent the approximately constant population of the pine 
trees vulnerable to infection; pi will be taken as the population of the southern 
pine beetle; and p2 that of its most important natural predator, the clerid beetle 
Thanasimus duhius. 

The total entropy brought into the ecosystem or carried out of it through 
one-body transport processes can be written as (Michaelian, 2005), 

a=l 

where the sum is over all n = 3 species and Pa is the population of species a. 

represents the average rate of exchange, or flow, of entropy with the external 
environment per individual per unit time of species a, as a result of energy 
(including heat), work, or matter flow (see below). 

Similarly, the internal entropy production, including individual production 
and exchange of entropy between individuals of the species, may be written as 



diS \ ^ 



dt 

a=l 



+ 0(4) 



a'=l 



> 0. 



(9) 

The Ta represent the entropy production of species a due to one-body irre- 
versible processes occurring within the individual such as; photosynthesis, tran- 
spiration, respiration, metabolism, etc. The r^Q:' represent the entropy pro- 
duction and exchange due to 2-body interactions between individuals of species 
a and a (e.g.. those involved in competition, predator-prey, symbiosis, etc.); 
Taa'a" Correspond to similar but 3-body interactions, and 0(4) represents the 
entropy production due to 4-body and higher order interactions (for example, 
those required for the functioning of societies). The 4-body and higher order 
n-body terms will be neglected in what follows since they would normally be 
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small as they require increasingly improbable (except for social species) n-body 
localization in space and time. 

Equation ([7]), the general evolutionary criterium, for the time change in the 
entropy production due to a change in the generalized forces X (the populations) 
then becomes 



dxV = ^ dpa 



Oi' a." 



< 0. (10) 



The dynamics of the ecosystem can now be determined from equations 
©, ® and Uni) once the T are specified. 

The r's represent entropy production and flow between individuals of the 
species and between individuals and their environment. A general expression 
for this entropy flow comes from the Gibb's equation and results from the flow 
of energy, work, and matter (Prigogine 1967). For example, the energy per 
individual per unit time taken in through photosynthesis dca, or the heat dqa 
per individual per unit time transported to the external environment, the work 
done on the environment per unit time PdVa at constant pressure P, and the 
matter components (e.g.. nutrients) of type /? taken in or given out by species 
a, duafs, of chemical potential /i^^, give for the rate of entropy exchange per 
individual with the environment, 

_ 1 deg + dqa P dVa _ ^ ST^ dugp 
T dt ^ T dt t2^^^P dt ' 

where the temperature T (of the participating individuals) may be approximated 
as being constant for the ecosystem (Gallucci, 1973). A similar expression can 
be written for the Taa' , representing the entropy production and exchange be- 
tween individuals of species a and a', i.e. in terms of the energy, work, and 
matter exchanged due to the 2-body interactions between individuals of the 
species. The affect on the entropy flows due to a simultaneous interaction of a 
third individual (three-body effects) of species a" is considered in the parameter 

Taa'a" ■ 

Determining the F's for a real ecosystem therefore requires the determination 
of the flows of energy c?e, heat dg, volume dV , and mass drip of type [3 between 
individuals of the participating species and between individuals and the exter- 
nal environment. Such a determination is possible in principle but obviously 
difficult in practice. Some of the experimental details for obtaining these types 
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of flow measurements between the individuals, and between the individuals and 
their environment, can be found in (Gallucci, 1973) and references therein. We 
have outlined the empirical determination of the entropy production coefficient 
Ta for a plant (Hernandez, 2008, Hernandez and Michaelian, 2008). The en- 
tropy exchange between individual animals Taa' is more difficult to measure. 
However, since much of the entropy flow between individuals is directly related 
to their diet, it is plausible, for example, that values for the Taa' may one day be 
obtained through a chemical analysis of the DNA of the excrement of individu- 
als. Similar analyses have in fact been performed for real ecosystems, however, 
with the focus on energy flow (see for example. Homer et al., 1976). 

In the absence of real ecosystem data concerning the production of entropy 
and the exchange of entropy, here we generate these coefficients F for a model 
ecosystem subject to thermodynamic law with the aid of a genetic algorithm 
(Michaelian, 1998). The algorithm begins by randomly generating sets of initial 
values for the F's within fixed ranges. The algorithm then evaluates the fitness 
of each set by checking to see if the population dynamics, as determined by 
criteria ([9]) and (fTO|) . leads to a viable stationary state in the long time limit, 
i.e. one with limt^oo diS/dt = —deS/dt (as required by classical irreversible 
thermodynamics for constant external constraints, see Eq. ([4])) and with diS/dt 
large and positive (consistent with what is known about the natural evolution of 
biotic systems to ever higher entropy production regimes (see Prigogine, 1967)). 
The best sets of F's are selected and evolved through mutation and crossover, 
optimizing (maximizing), the fitness function, 

d,s/dt ^^2^ 



d^S/dt + deS/dt' 

which, as required, is large for dtS/dt large and for diS/dt = —deS/dt. 

An example of a set of F's so obtained is given in Appendix A. Using this set 
of F's and a starting population set of populations {pi — 1000, P2 = 2000) with 
P3 fixed at 2000 (the constant external constraint), and generating infinitesimal 
variations of the populations dpi and dp2 at random {dp^ = 0) while only 
accepting those sets {dp} which satisfy the thermodynamic criteria of Eqs. ([9]) 
and (fTO| . leads to the stable cyclic attractor (oscillating) population dynamics 
as shown in figure [TJ 
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Figure 1: (a) Populations pi (solid line) and p2 (dashed line) as a function of 
time k for the ecosystem given in Appendix A in response to various (see text) 
perturbations, (b) diS/dt (solid Hne) and —deS/dt (dashed line) as a function of 
time, (c) Trajectory in population space showing the cychc attractor dynamics 
and the effects of 4 distinct perturbations (see text). 
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IV. Response to Perturbation 

The response to perturbation obtained under the dictates of the thermody- 
namic criteria, Eqs.. ([9]) and (fTO|) . once the ecosystem has arrived at a stationary 
state, Eq. ([5]), are shown in figure [TJ This stationary state, for the particular 
set of r's obtained, is a cychc attractor. The first perturbation, affected at 
time k = 5e6 by reducing the population of pi to 40, can be seen to have little 
effect on the ecosystem, a full recovery of the populations is obtained rapidly. 
The second perturbation affected by increasing the population of pi to 800 
at time k = 8e6 also produces only a small transient effect. However, if the 
ecosystem is perturbed in the same manner but at time k = 1.3e7, population 
Pi goes negative (as does the internal production of entropy diS/dt) implying 
the extinction of the species and the non-viability, thermodynamically, of the 
ecosystem. The timing of a perturbation affected on an ecosystem with cyclic 
attractor population dynamics therefore appears to be crucial in deciding the 
fate of the ecosystem. 

Similar results are obtained if population p2 is perturbed. A perturbation 
of the ecosystem at time k = 1.7e7 by increasing population p2 to 1600 has no 
long term consequences. However, when population p2 is reduced to 600 at time 
k = 2.0e7, the internal production of entropy and the population pi both go 
negative, the ecosystem again becomes untenable. 

Figure [2] plots the dynamics of the ecosystem in population space pi : p2 for 
50 different initial populations. It is apparent from this figure that the recovery 
or not of an ecosystem from a particular perturbation depends on the region in 
population space into which the ecosystem is perturbed. Perturbation into the 
"regions of danger" marked on the figure leads to either one of the populations 
going negative, or to the internal production of entropy going negative. Both of 
these results are non physical and would foretell the collapse of the ecosystem 
in nature. Interestingly, these regions of danger do not necessarily correspond 
to regions of small population. 

We next examine the response of our model ecosystem to permanent changes 
in the boundary conditions. Figure [3] shows the ecosystem dynamics under new 
boundary conditions olpz reduced to 200, down from its originally fixed value of 
P3 = 2000, implying a less negative fiow of entropy into the ecosystem. Without 
allowing time for the interaction coefficients F to evolve in response to the new 



11 



Figure 2: The dynamics in population space for 50 different initial populations 
showing that there are "regions of danger" in population space for which a 
perturbation into these regions would cause the system to become untenable, 
either because one of the populations extinguishes (goes negative) or because 
the internal production of entropy becomes negative (non-physical). 
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boundary conditions, the cyclic attractor then becomes a point attractor as 
shown in figure [3l However, as can be seen from this figure, the point attractor 
is not a thermodynamic stationary state since the internal production of entropy 
is no longer equal to the negative of the external flow of entropy (equation ^ 
is no longer satisfied). The ecosystem fitness function, Eq.. (|12p. is no longer 
at a local maximum value and the system, given time, would evolve its entropy 
production and exchange coefficients (the set F) until reaching a new stationary 
state where the production and external flow of entropy are once again equal. 
Note that in the perturbed state with the new external constraint, = 200, the 
same small perturbation of reducing the population pi to 40 at time k = 5e6, 
which did not have any lasting affect on the ecosystem previously (Fig. [1]), now 
results in the collapse of the ecosystem since it moves it into the non-physical, 
thermodynamically prohibited, regime of negative internal entropy production 
(figure [S]). 

Figure |4] shows the opposite effect of increasing the flow of negative entropy 
into the ecosystem, obtained by increasing the value of the fixed external con- 
dition to p3 = 2200. Without allowing time for the interaction coefficients F 
to evolve, the dynamics remains that of a cyclic attractor but the orbit of the 
attractor increases significantly, bringing the population pi very close to zero at 
one point in its orbit. A slight perturbation of population p2 at time k = 5e6 
is sufficient to cause the population pi to pass through zero and thereby cause 
the collapse of the ecosystem. We believe that this is a possible thermodynamic 
explanation of the "enrichment paradox" (Rosenzweig, 1971); contrary to naive 
expectation, an increase in the inflow of nutrients is often observed to make an 
ecosystem more vulnerable to perturbation. This will be considered in detail in 
a forthcoming paper (Alonso, 2007). 

We have also verified (Alonso, 2007) that had the set of F's been chosen such 
that the thermodynamic stationary state corresponded to a point attractor in 
population space (sometimes referred to as an equilibrium state in the ecological 
literature), then increasing the value of the fixed external condition p^ leads to 
population oscillations of pi and p2. This may have relevance to the sudden 
outbreak of large population oscillations of certain insect pathogens, such as 
the southern pine beetle, as described in the following section. 
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Figure 3: Ecosystem dynamics for the case where the fixed boundary conditions 
have been changed from ps = 2000 to = 200. The population dynamics be- 
comes that of a point attractor. However, the system is not in a thermodynamic 

stationary state since diS/dt ^ —deS/dt. A subsequent perturbation of pi to 
40 at fc = 5e6, which did not affect the stabihty previously (see figure 1), now 
leads to diS/dt going negative, violating the 2nd law. The vulnerability of the 
ecosystem has thus been increased by reducing the external constraint of ps . 
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Figure 4: Ecosystem dynamics for the case where the fixed boundary condition 
has been increased from ps = 2000 to ps = 2200. In this case, the population 
dynamics is that of a cychc attractor with a larger orbit, bringing pivery close 
to zero at times. A subsequent small perturbation of p2 to 800 at fc = 5e6 leads 
to diS/dt going negative and population pipassing through zero. We believe 
that this is the thermodynamic origin of the "enrichment paradox" (Rosenzweig 
1971). 
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V. Southern Pine Beetle Population Oscillations 



The Southern Pine Beetle, indigenous to the southern and south eastern 
states of the United States and to the northern parts of Mexico, is one of the 
pathogens of most concern to humans as it can severely affect the economic value 
of a forest stand. Its population dynamics often exhibits infestation outbreaks 
of an oscillatory nature. This dynamics has been extensively studied with a 
view to understanding its origin and thereby to controlling the pest (Turchin et 
al., 1999). In the wild state of natural forests, the pine beetle populations are 
small and relatively constant and the insect may be considered as an integral 
component of a robust ecosystem, even helping to promote normal ecosystem 
succession (Natural Resources Canada, 2003). However, the last century has 
seen an almost complete replacement of the natural mixed forest stands with 
homogeneous stands of particular species of greater economic value to humans, 
such as soft wood pine. This has led to forests more vulnerable to pine beetle 
attack with the result that millions of hectors of pine trees can be devastated by 
a single beetle outbreak (Aukema et al., 2008). The outbreak is often preceded 
by a debilitating event affecting the forest such as a drought, fire, or flood 
(McNulty et al., 1998). 

The details surrounding an outbreak of a sustained pine beetle attack are 
well known. Normally the pine beetle population is relatively small and constant 
with the beetle preying only on the vulnerable weaker trees of a stand. The 
vulnerable pine trees provide nourishment and protection for the pine beetle 
which buries under the inner bark of the tree to deposit its eggs in bored out 
gallerias. The beetles carry a fungus which quickly establishes itself within 
the tree and spreads out to cut off the water and nutrient irrigation of the 
tree. The vulnerable tree normally dies within a year of first attack. Healthy 
trees efficiently repel the pine beetle attack by emitting enough resin in the 
beetle galleries to dissuade further colonization (Blinka, 2007) and induce beetle 
mortality through a number of effects including crystallization of the sap around 
beetle larva. 

Unlike vulnerable pine trees, healthy pine trees, therefore, cannot normally 
be considered as contributing to the primary host species of the pine beetle. If, 
however, there occurs a debilitating event, such that the density of vulnerable 
pine trees in a given area becomes great enough, then the local population of 
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pine beetles may swell to such a point that there is now a large enough number 
of beetles on each tree 200/m^) that even the healthy trees cannot cope and 
become vulnerable (Blinka, 2007). 

It has long been debated as to whether exogenous (originating externally) or 
endogenous (originating from within) factors contribute most to the observed 
population oscillations of the pine beetle. Turchin (Turchin, 2003) has argued for 
a synthetic approach, suggesting that both enogenous and endogenous factors 
can be important. For example, while it is known that particularly cold winters 
can reduce and even end an outbreak, otherwise very little correlation of climate 
on population dynamics has been found (Turchin ct al., 1991). This fact has 
lead to a search for possible endogenous effects to explain oscillations in the 
population density of the pine beetle. 

One endogenous effect studied by Turchin et al. is a time delay in density- 
dependent affect of a predator on the pine beetle. The effect of introducing 
a time delay into the normally stabilizing density dependent terms of Lotka- 
Volterra type equations was first studied theoretically by May (1974) for the 
systems of vegetation-herbivore, and vegetation-herbivore-predator. For the 
vegetation-herbivore system, May found that a time delay in the normally sta- 
bilizing density term of the logistic equation could lead to limit cycles (stable 
population oscillations) if the characteristic growth rate of the herbivore was 
small compared to the time delay in the regulatory mechanism (stabilizing term 
of the logistic equation related to the density of the herbivore). When a preda- 
tor was included to prey on the herbivore. May found it much more difficult 
to obtain oscillatory solutions to the model equations. Nevertheless, through 
carefully controlled field experiments, Turchin et al. (1999) verified that there 
indeed exists a delayed density-dependent affect of predators on the pine beetle 
which could thus contribute to the population oscillations. 

The essence of the pine beetle ecosystem can be reduced to three strongly 
interacting species; the primary production species (the vulnerable pine tree), 
the herbivore (the pine beetle itself) and a predator species (the clerid beetle, 
the most important natural predator of the pine beetle). This reduction of the 
ecosystem to only three species is convenient (but not necessary) to demonstrate 
how our formalism, based on the application of general thermodynamic criteria 
concerning entropy production and exchange, can describe the dynamics of the 
pine beetle ecosystem. 
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In section HVl we have shown that if the interactions (exchange and produc- 
tion of entropy) between individuals of the different species of an ecosystem is 
of n-body nature (including higher than 2-body terms) then a stable point at- 
tractor or stable limit cycle dynamics can be obtained after the system relaxes 
to its thermodynamic stationary state (see figure ??). Consider now such a 
stable ecosystem in its stationary state, being either a point attractor or stable 
limit cycle in population space. The initial constant boundary conditions for 
this ecosystem correspond to that segment of the total population of pine 
trees within a given area that are vulnerable to attack while the ecosystem is 
enjoying normal, stable conditions. If the population of the vulnerable trees 
then increases suddenly due to a debilitating event affecting the forest, then 
the beetle population pi and the Clerid beetle predator population p2 start to 
oscillate with potentially much greater amplitude (Fig. 2]). 

The interaction (production and exchange of entropy) of the pine beetles 
with the pine trees and its predator is non-linear in their respective populations 
and the physical restrictions on the total entropy production (the 2nd law) and 
on the sign of the rate of change of the production of entropy (Prigogine's gen- 
eral evolution criterium) together force this oscillatory population dynamics. 
However, this is a perturbed ecosystem, one which is no longer in a thermody- 
namic stationary state, as could be verified in the field by demonstrating that 
the production of positive entropy has temporarily increased and is no longer 
equal to that of the inflow of negative entropy (see Eq. 2] and figured]). 

This explanation of the population oscillations is in accord with the follow- 
ing known facts regarding pine beetle population dynamics; 1) the pine beetle 
populations are not always oscillating but can be relatively constant in healthy 
forest stands, 2) oscillations usually begin after a particular debilitating event 
which causes a rapid increase in the number of vulnerable trees, the primary 
resource available to the beetle, 3) the amplitude of the population oscillations 
are correlated with the increase in the amount of vulnerable trees in a given 
area, or, in other words, to the density of nominally healthy pine trees in the 
area, 4) it is known that perturbed ecosystems (including diseased systems) have 
temporarily increased total entropy production (Schneider and Sagan, 2005), 5) 
it is contingent on both exogenous (debilitating events) and endogenous (ther- 
modynamic laws) factors. 

The results presented above suggest a possible means for controlling a pine 
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beetle outbreak. Figure [5] shows that there are regions in population space for 
which perturbation into these regions leads to eventual collapse of the ecosystem. 
Once the entropy flow coefficients T have been determined, the thermodynamic 
formalism presented here determines these regions of danger. It is then only a 
matter of perturbing one of the populations into these regions of danger to obtain 
the desired extinction of the outbreak. This could be obtained by augmenting, 
or reducing, the population of one of the involved species at the indicated point 
in the population cycle. 

VI. Summary and Conclusions 

Acknowledging that ecosystems, like all macroscopic processes, are subject 
to definite thermodynamic law, we have demonstrated that under constant ex- 
ternal constraints, the response of ecosystems to perturbation can, in principle, 
be predicted. The thermodynamic criteria which direct the dynamics come 
from non-equilibrium thermodynamic theory. They are; 1) the system must 
eventually arrive at a thermodynamic stationary state, Eq. 2) the internal 
production of entropy must be positive definite, in accord with the second law 
of thermodynamics, Eq. and 3) any natural change in the populations must 
be in such a manner so as to reduce the internal production of entropy of the 
entire system, Prigogine's general evolution criterion, Eq. ([7|). 

In the absence of data on the entropy production and exchange of real ecosys- 
tems, we considered a simple model ecosystem generated by evolving interaction 
coefficients (representing entropy production and exchange) through selection 
with a fitness function favoring the thermodynamic criteria identified above. 
We then studied the response of this model ecosystem to perturbation of the 
populations under the same thermodynamic criteria. We found that there exists 
regions in population space for which perturbation into these regions leads to 
the eventual extinction of one or more of the species, or to a negative internal 
production of entropy. The latter violates the second law of thermodynamics 
and would lead to ecosystem collapse since physical maintenance processes re- 
quire positive production of entropy. An important finding is that these regions 
in population space are a general feature of the thermodynamic framework and 
do not necessarily correspond to regions of small population. Assigning species 
to "in danger of extinction" lists, solely on the basis of the smallness of their 
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populations may therefore not be an effective conservation policy. Our proposed 
approach, based on thermodynamic criteria, predicts the dynamics over all of 
population space and thus leads to quantitative elements for providing more 
informed policy statements for responding to ecosystem perturbation. 

Increasing or decreasing the negative flow of entropy (natural resources) into 
the ecosystem has the effect of increasing or decreasing respectively both the 
amplitude of the orbit of the attractor in population space and the internal 
production of entropy of the system. In either case, this results in a more 
vulnerable ecosystem since the populations pass closer to zero or the internal 
production of entropy may more easily become negative respectively. We believe 
that this result is a thermodynamic explanation of the "enrichment paradox" . 

We applied our thermodynamic framework to a reduced ecosystem consisting 
of the southern pine beetle, pine trees, and the beetles most important predator. 
We have shown that the population oscillation of the southern pine beetle may 
be viewed within this thermodynamic framework as resulting from an increase 
in the inflow of resources (negative entropy) into a system which has a non-linear 
dependence of the production and flow of entropy on the species populations. 
Our thermodynamic description is contingent on the initial conditions known 
to precede an attack (the increase in the density of the vulnerable trees), pre- 
dicts pine beetle population oscillations for perturbed ecosystems, and predicts 
constant (or small amplitude oscillations) populations for unperturbed ecosys- 
tems. The model is consistent with the synthetic view that oscillations of the 
pine beetle populations are a result of both exogenous factors (the initial debil- 
itating event causing an increase in the number of vulnerable pine trees) and 
endogenous factors (thermodynamic laws and non-linear coupling of entropy 
flows). 

We proposed a method of control of the pine beetle by first delimiting these 
regions of danger, and then perturbing the populations into these regions. 

Ecosystems are composed of many thousands of interacting species and the 
details of the dynamics is, undoubtedly, signiflcantly more complicated. How- 
ever, our thermodynamic framework can be straightforwardly applied to a much 
larger and more complex ecosystem simply by measuring all the entropy pro- 
duction and exchange coefficients F for all the species involved. Work in this 
direction is underway (Hernandez, 2008). 

In conclusion, as for all macroscopic process, ecosystems are subject to def- 
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inite thermodynamic law. For constant external constraints, these laws are 
sufficient to determine ecosystem response to perturbation. Our analysis of the 
population dynamics based on thermodynamic law and the formulation of the 
interaction coefficients in terms of physical and measurable quantities (the pro- 
duction and exchange of entropy) is a step toward a more quantitative theory 
of ecosystems. 
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Appendix A 



The following set of species interaction (entropy production and exchange) co- 
efhcients (see Eqs. (|5]) and for an n = 3 species ecosystem and including 
up to 3-body terms, were used for calculating the dynamics of figures [T] to SJ 



= 1.1721 
Ti = 0.6277 
Til = -0.01616 

= 0.0372 

= 0.01617 
Till = 0.000448 
ri22 = -0.0000574 
rii2 = -0.0001446 



r| = 0.4710 
r2 = 8.787 
T22 = -0.04886 
r2i 0.0222 
r23 = 0.00 3 39 
r222 = -0.0000078 
T211 = -0.0002856 
r22i = 0.0002881 



r| = -55.000 
Tg = 1.1884 
Tgs = 0.020869 
r32 = 0.012619 



These coefficients were obtained by evolving various initial random sets of co- 
efficients, selected from within a finite range, using a genetic algorithm with a 
fitness function favoring a stationary state, i.e. diS/dt = —d^S/dt and diS/dt 
large in the long-time limit, (ie. at the stationary state). Although our 3 species 
ecosystem is artificial, the coefficients have characteristics of what would be ex- 
pected of a real ecosystem. For example, the coefficient of entropy fiow with the 
external environment, r|, representing the fixed external constraints, is of large 
negative value. The one-body internal entropy production coefficients Fi . . . Fa 
are all positive (respecting the 2nd law of thermodynamics) and of large value 
with respect to the 2-, and 3-body coefficients. The second law for the to- 
tal internal entropy production, Eq.([2|), is also respected over a large region of 
population space (see figure [2]). 
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